active_sites = [1 2 4 7 8];
passive_sites = [3 5 6 9];
control_sites = 10;

for cur_exp = 1:9
    for layer = 1:1
        for time_point = 1:4
            cur_fb_act = proj_meta(cur_exp).rd(layer,time_point).act(:,1:7500);
            cur_pb_act = proj_meta(cur_exp).rd(layer,time_point).act(:,7501:15000);
            cur_fb_vel = -proj_meta(cur_exp).rd(layer,time_point).velM_smoothed(1:7500);
%             cur_fb_P = proj_meta(cur_exp).rd(layer,time_point).velP_smoothed(1:7500);
            cur_pb_vel = -proj_meta(cur_exp).rd(layer,time_point).velM_smoothed(7501:15000);
%             cur_pb_P = proj_meta(cur_exp).rd(layer,time_point).velP_smoothed(7501:15000);

            % generate binary activity during feedback.
            velM_fb_bi = proj_meta(cur_exp).rd(layer,time_point).velM_smoothed(1:7500);
            velM_fb_bi = -velM_fb_bi; % make the velocity of forward-running positive.
            velM_fb_bi(velM_fb_bi > 0.001) = 1;
            velM_fb_bi(velM_fb_bi <= 0.001) = 0;
            
            % generate binary activity during playback.
            velM_pb_bi = proj_meta(cur_exp).rd(layer,time_point).velM_smoothed(7501:15000);
            velM_pb_bi = -velM_pb_bi; % make the velocity of forward-running positive.
            velM_pb_bi(velM_pb_bi > 0.001) = 1;
            velM_pb_bi(velM_pb_bi <= 0.001) = 0;
            
            % generate the frames of ps onset.
            ps_ons = find(diff(proj_meta(cur_exp).rd(layer,time_point).ps_id(1:7500))>1);
            % remove early & late ps.
            ps_ons(ps_ons > 7425) = [];
            ps_ons(ps_ons < 6) = [];
            
            % generate the frames of playback halt onset.
            pb_ons = find(diff(proj_meta(cur_exp).rd(layer,time_point).ps_id(7501:15000))>1);
            % remove early & late pb.
            pb_ons(pb_ons > 7425) = [];
            pb_ons(pb_ons < 8) = [];
            
            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(ps_ons)
                if sum(velM_fb_bi(ps_ons(ps_id)-5:ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(cur_exp,time_point) = rps_nbr;
                end
            end
            
            rpb_nbr = 0; % number of real pb halt
            rpb_id = [];
            for pb_id = rps_id
                if sum(velM_pb_bi(pb_ons(pb_id)-7:pb_ons(pb_id)+22)) == 0
                    rpb_nbr = rpb_nbr+1;
                    rpb_id(rpb_nbr) = pb_id;
                    all_rpb_nbr(cur_exp,time_point) = rpb_nbr;
                end
            end
            
%             if rps_nbr > 10
%                 rps_id = rps_id(1:10);
%             end
%             
%             if rpb_nbr > 10
%                 rpb_id = rpb_id(1:10);
%             end
            
%             for ROI_id = 1:size(cur_fb_act,1) % mm_cells{cur_exp,layer} % change here if calculation of selected cells needed
            rps_resp_cnt = 0;
            rps_resp = [];
            for rps_act_nbr = 1:length(rps_id)
                rps_frame_id = ps_ons(rps_id(rps_act_nbr));
                rps_resp_cnt = rps_resp_cnt+1;
                rps_resp(rps_resp_cnt,:) = cur_fb_vel(rps_frame_id-7:rps_frame_id+75);
            end
            rps_resp_by_cell{cur_exp,time_point} = rps_resp;
%             end
            
%             for ROI_id = 1:size(cur_pb_act,1) %mm_cells{cur_exp,layer} % change here if calculation of selected cells needed
            rpb_resp_cnt = 0;
            rpb_resp = [];
            for rpb_act_nbr = 1:length(rpb_id)
                rpb_frame_id = pb_ons(rpb_id(rpb_act_nbr));
                rpb_resp_cnt = rpb_resp_cnt+1;
                rpb_resp(rpb_resp_cnt,:) = cur_pb_vel(rpb_frame_id-7:rpb_frame_id+75);
            end
            rpb_resp_by_cell{cur_exp,time_point} = rpb_resp;
%             end
        end
        % comment out this line if calculation of selected cells needed
%         cell_cnt{cur_exp,layer} = size(proj_meta(cur_exp).rd(layer,1).act,1);
    end
end

% cell_cnt = mm_cells_cnt; % and here
for cur_exp = 1:9
    for layer = 1:1
%         if cell_cnt{cur_exp,layer} > 0
%         for ROI_id = 1:cell_cnt{cur_exp,layer} % mm_cells{cur_exp,layer} % and here
        cur_fb_cell_act = vertcat(rps_resp_by_cell{cur_exp,:});
        if size(cur_fb_cell_act,1) == 1
%             cur_fb_cell_act = cur_fb_cell_act-mean(cur_fb_cell_act(1:7));
            mean_rps_resp_by_cell{cur_exp} = cur_fb_cell_act;
        else
            mean_cur_fb_cell_act = mean(cur_fb_cell_act);
%             mean_cur_fb_cell_act = mean_cur_fb_cell_act-mean(mean_cur_fb_cell_act(1:7));
            mean_rps_resp_by_cell{cur_exp} = mean_cur_fb_cell_act;
        end
%         end
%         end
    end
end

for cur_exp = 1:9
    for layer = 1:1
%         if cell_cnt{cur_exp,layer} > 0
%             for ROI_id = 1:cell_cnt{cur_exp,layer} %mm_cells{cur_exp,layer} % and here
        cur_pb_cell_act = vertcat(rpb_resp_by_cell{cur_exp,:});
        if size(cur_pb_cell_act,1) == 1
%             cur_pb_cell_act = cur_pb_cell_act-mean(cur_pb_cell_act(1:7));
            mean_rpb_resp_by_cell{cur_exp} = cur_pb_cell_act;
        else
            mean_cur_pb_cell_act = mean(cur_pb_cell_act);
%             mean_cur_pb_cell_act = mean_cur_pb_cell_act-mean(mean_cur_pb_cell_act(1:7));
            mean_rpb_resp_by_cell{cur_exp} = mean_cur_pb_cell_act;
        end
%             end
%         end
    end
end

resp_fb_act = vertcat(mean_rps_resp_by_cell{active_sites});
resp_fb_pas = vertcat(mean_rps_resp_by_cell{passive_sites});

mean_resp_fb_act = mean(resp_fb_act);
mean_resp_fb_pas = mean(resp_fb_pas);

sem_resp_fb_act = std(resp_fb_act)/sqrt(size(resp_fb_act,1));
sem_resp_fb_pas = std(resp_fb_pas)/sqrt(size(resp_fb_pas,1));

figure;hold on
h1 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_fb_act+sem_resp_fb_act fliplr(mean_resp_fb_act-sem_resp_fb_act)]);
h2 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_fb_pas+sem_resp_fb_pas fliplr(mean_resp_fb_pas-sem_resp_fb_pas)]);
plot(-7/15:1/15:75/15,mean_resp_fb_act)
plot(-7/15:1/15:75/15,mean_resp_fb_pas,'r')
xlim([-0.5 4.5])

resp_pb_act = vertcat(mean_rpb_resp_by_cell{active_sites});
resp_pb_pas = vertcat(mean_rpb_resp_by_cell{passive_sites});

mean_resp_pb_act = mean(resp_pb_act);
mean_resp_pb_pas = mean(resp_pb_pas);

sem_resp_pb_act = std(resp_pb_act)/sqrt(size(resp_pb_act,1));
sem_resp_pb_pas = std(resp_pb_pas)/sqrt(size(resp_pb_pas,1));

h3 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_pb_act+sem_resp_pb_act fliplr(mean_resp_pb_act-sem_resp_pb_act)]);
h4 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_pb_pas+sem_resp_pb_pas fliplr(mean_resp_pb_pas-sem_resp_pb_pas)]);
plot(-7/15:1/15:75/15,mean_resp_pb_act,'--')
plot(-7/15:1/15:75/15,mean_resp_pb_pas,'r--')
xlim([-0.5 4.5])
% ylim([-1 3.5])

set(h1,'FaceColor',[.5 .5 .5])
set(h2,'FaceColor',[.5 .5 .5])
set(h3,'FaceColor',[.5 .5 .5])
set(h4,'FaceColor',[.5 .5 .5])
set(h1,'LineStyle','none')
set(h2,'LineStyle','none')
set(h3,'LineStyle','none')
set(h4,'LineStyle','none')